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AEROSPACE 


A  Multifunction  Space-Borne  Scanner 

90FE0179A  Beijing  YUHANG  XUEBAO  [JOURNAL  OF  CHINESE  SOCIETY  OF  ASTRONAUTICS] 
in  Chinese  No  2,  30  Apr  90  pp  52-56 

[Article  by  Kong  Huixing  [7895  1920  5281],  Shanghai  Institute  of  Technical 
Physics] 

[Text]  Abstract 

A  proposal  is  made  to  add  five  detection  bands  to  an  existing  very-high- 
resolution  (VHR)  scanning  radiometer  to  enhance  its  capability  for  the 
observation  of  meteorological  phenomena,  ocean  environment,  and  natural 
disasters. 

(I)  Performance  of  the  Existing  Very-High-Resolution  Scanning  Radiometer 

Test  results  show  that  the  performance  indices  of  China's  existing  VHR 
scanning  radiometer  in  the  visible-light  and  infrared  band  meet  all  the 
design  requirements,  as  shown  in  Table  1.  Both  visible-light  and  infrared 
channels  produce  clear  images,  and  have  good  observation  capability  for 
cloud  maps,  ground  covers,  lakes  and  ocean  water. 

(II)  Satellite  and  Radiometer  Design  Requirements  for  Monitoring  the 
Environment  and  Natural  Disasters 

Because  natural  disasters  occur  abruptly  in  time  and  in  space,  real-time 
monitoring  is  required  in  order  to  provide  timely  reporting  and  assistance. 
An  effective  approach  of  monitoring  the  entire  Chinese  landmass  is  to  use  a 
satellite-bome  VHR  scanning  radiometer  operating  at  high  repetition 
frequency . 

The  results  of  recent  studies  in  earth  science  and  environmental  science 
show  that  in  order  to  explain  the  evolution  of  certain  important  environ¬ 
mental  phenomena,  it  is  necessary  to  consider  the  earth  as  a  single  entity, 
and  to  study  the  global  effects  between  the  atmosphere,  the  ocean,  and  the 
land,  the  gain  and  loss  of  earth  radiation,  and  the  combustion  of  fossil 
fuels,  etc.  To  accomplish  this  requires  the  accumulation  of  over  10  years 
of  observations  with  high  spatial  and  time  resolution,  as  indicated  in 
Table  2. 
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Table  1.  On—Orbit  Performance  of  Very— High— Resolution  Scanning  Radiometer 

Parameter 

Design  specification  j  Test  results 

Detection 

band 

i 

t ''  . . .  .  . . .  ■  . 

0.58-0.68  ym  1  Visible  cloud  maps,  ground  cover,  cloud 

I  and  snow  separation,  ocean  sediment 
0.725-1.1  ym  i  Visible  cloud  maps,  ground  cover,  cloud 

j  and  snow  separation,  water  and  land 

* 

1  0.48-0.53  ym 

0.53-0.58  ym 

j  10.5-12.5  ym 

i 

border  separation 

Ocean  water  color,  chlorophyll 
distribution 

Ocean  water  color,  chlorophyll 
j  distribution 

;  Infrared  cloud  maps,  water  and  land 
border  separation,  temperature 

I  measurement 

Spatial 

resolution 

j  1.08  km 

1.08  km 

Detection 

sensitivity 

Ch^,  Ch2,  NEAp<0.16% 
Ch3,  Ch4,  NEAp<0.25% 
Ch5  NEAT£0.8  K 

^  NEAp<0.11% 

NEAp<0.15% 

NEAT<p.5  K 

Detection 

accuracy 

Ch3,Ch4  ±10% 

Ch5  1.0  K 

±10% 

approximately  1.0  K 

Table  2.  Time  and  Spatial  Resolution  Requirements  for  Environmental 
Satellite  Observations 


Obj  ect  of 
observation 

Repetition  period 

Ground 

resolution 

Atmosphere 

At  least  0.5  day /observation 

1  km 

Ocean 

0.5-1  day/observation  (for  ocean  conditions) 

3-5  days/observation  (for  medium  and  large-scale 

0.5-1  km 

ocean  phenomena) 

1-4  km 

Ground  cover 

7-10  days/observation 

0.5-1  km 

Although  earth  resource  satellites  have  very  high  spatial  resolution,  it 
cannot  satisfy  the  repetition-frequency  requirement  for  environmental 
monitoring  because  it  takes  16-26  days  to  cover  the  whole  earth.  In  addi¬ 
tion,  the  multi-band  scanner  of  an  environmental  satellite  has  5-10  times 
higher  sensitivity  than  that  of  an  earth  resource  satellite  (see  Table  3)  and 
therefore  can  detect  minute  changes  in  the  environmental  parameters  and 
predict  its  future  development .  Because  of  these  unique  features  of  the 
©nvironmental  satellite,  it  should  be  further  developed  to  provide  the 
capability  for  monitoring  the  environment  and  natural  disasters. 
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Table  3.  Comparison  of  Multi-Band  Scanners  of  an  Environmental  Satellite  and  of  an  Earth  Resource 
Satellite 
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(Ill)  10-Band  Scanning  Radiometer 

In  order  to  enhance  the  monitoring  capability  of  the  existing  VHR  scanning 
radiometer,  it  is  necessary  to  add  five  more  detection  bands  to  the  existing 
bands  to  provide  a  10-band  observation  capability.  Specifically,  the 
enhancements  include  the  following: 

(1)  Adding  a  3.55-3.93  ym  medium-infrared  band  for  monitoring  forest  fires 
and  grassland  fires;  the  low  water  absorption  in  this  band  will  provide 
improved  temperature  measurement  accuracy  over  the  ocean  surface  in  the  mid- 
low  latitude  region. 

This  band  has  very  high  sensitivity  with  respect  to  fire  because  it  is 
located  at  the  radiation  peak  of  a  780  K  target.  The  following  table  shows 
the  calculated  results  of  radiation  luminance  for  targets  at  different 
temperatures . 


Table  4.  Radiation  Luminance  of  Targets  at  Different  Temperatures  in  the 
10.5—12.5  ym  Band  and  the  3.55—3.93  ym  Band 


la^et  temperature 

T=1100  K 

T=800  K 

! 

'  T=300  K 

I  N(1100  K)  1 

N(800  K) 

Target  liaminance^^'''-'--^ 

N(300  K) 

N(300  K) 

3.55-3.93  ym 
Luminance  (W/cm2 • ster) 

0.1924 

5.08x10-2 

1.7x10"5 

1 

1.13x10^ 

3.9x103 

10.5-12.5  ym 
Ltiminance  (W/ cm2 .  ster) 

5.59x10-2 

3.13x10-2 

1.86x10-3 

30.0 

16.8 

The  radiated  power  received  by  a  scanner  from  a  target  which  does  not  fill 
the  instantaneous  field  of  view  is  proportional  to  the  target  area.  Thus, 
the  radiated  power  of  a  20x20  m^  fire  at  800®K  is  equivalent  to  that  of  a 
!• 1x1.1  Km  target  at  300°K.  Taking  into  account  atmospheric  attenuation 
and  margin  threshold  setting,  it  is  estimated  that  the  scanner  can  detect  a 
30x30  m2  open  fire  and  issue  a  fire  alarm. 

(2)  Adding  a  1.58-1.64  ym  shortwave  infrared  band  for  measuring  the  water 
content  in  soil  and  crops  and  for  monitoring  drought  conditions. 

The  soil  reflectivity  drops  sharply  in  the  1.3-2. 5  ym  band  when  it  contains 
water;  in  particular,  there  are  water  absorption  peaks  at  1.44  ym,  1.90  ym 
and  2.2  ym.  When  the  soil  moisture  is  in  the  0.8  percent  to  20.2  percent 
range,  an  increase  of  1  percent  in  water  content  will  reduce  the  reflec¬ 
tivity  by  1  percent  (see  Figure  1) .  The  water  content  in  the  leaves  of 
ground  vegetation  will  clearly  affect  the  reflectivity  in  the  1.3-2. 5  ym 
band  (see  Figure  2) .  The  sensitivity  of  the  scanning  radiometer  is  estimated 
to  NEAp£0.1%,  which  is  sufficient  to  detect  a  1  percent  change  in  water 
content;  therefore,  the  scanner  is  effective  in  monitoring  drought  condi¬ 
tions.  Considering  the  fact  that  there  are  H2O  and  CO2  absorption  in  the 
atmosphere,  the  detection  band  is  chosen  to  be  1.58-1.64  ym. 
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Figure  1.  Spectral  Characteristics  of  Soil  Reflectivity  With 
Different  Moisture  Content 


Figure  2.  Reflectivity  Spectrum  of  Clay  and  Wheat  Leaves  With 
Different  Moisture  Content 

(3)  Adding  a  0.433-0.453  ym  band  for  measuring  the  chlorophyll  content  in 
the  ocean.  This  band  supplements  the  existing  0.48-0.53  yin  band  and  the 
0.53-0.58  ym  band  to  measure  the  water  color  which  reflects  a  wide  range  of 
chlorophyll  content.  The  collected  data  can  not  only  be  used  for  preliminary 
survey  of  ocean  productivity  but  also  be  used  in  conjunction  with  surface 
temperature  data  for  fishing-ground  forecast. 

(4)  Because  of  scattering  by  molecules  and  colloidal  particles  in  the 
atmosphere  in  the  0.433-0.453  ym  band,  the  sky  background  is  one  order  of 
magnitude  brighter  than  the  ocean  surface.  In  order  to  Isolate  the  water- 
color  signal  of  the  ocean  from  the  signal  of  the  bright  sky  background,  it 
is  necessary  to  add  a  0.745-0.785  ym  band  to  correct  for  the  effect  of  the 

colloidal  particles  in  the  atmosphere. 

(5)  The  existing  10.5-12.5  ym  thermal  infrared  band  of  the  scanning 
radiometer  is  split  into  two  bands:  10.5-11.5  ym  and  11.5-12.5  ym,  in  order 
to  partially  cancel  the  atmospheric  effect  so  that  higher  measurement 
accuracy  of  the  ocean  temperature  (approximately  1.5  K)  can  be  achieved. 

(6)  Because  of  the  large  differences  between  the  reflectivities  of  soil, 
ground  cover,  and  water  in  the  0.725-1.1  ym  band,  it  is  an  important  band 
for  monitoring  land-water  separation,  river  flow  rate,  and  flood  damage. 
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In  order  to  improve  the  capability  of  separating  water  and  land,  the  ground 
resolution  in  this  band  is  increased  to  540  m. 

(7)  Currently,  the  0.725-1.1  ym  band  and  the  0.58-0.68  ym  band  are  combined 
to  monitor  the  ground  cover  index,  which  is  a  measure  of  ground  cover 
distribution,  growth  trend  of  crops,  and  estimated  wheat  production.  In 
order  to  match  the  spatial  resolution  between  these  two  bands,  the  resolution 
in  the  0.58-0.68  ym  band  is  also  increased  to  540  m. 


Figure  3.  Upward  Spectral  Luminance  of  Ocean  Water  With 
Different  Chlorophyll  Content 


I  Ref lectivity  % 
50- 


Figure  4.  Reflectivity  Spectrum  of  Typical  Targets 

The  modified  scanning  radiometer  has  10  detection  bands.  It  includes  the 
existing  five  bands  of  the  AVHRR  of  the  NOAA  weather  satellite  and  the 
1.58-1.64  ym  band  which  is  to  be  added  on  the  NOAA-K,  L,  M  satellites. 

The  addition  of  the  0.433-0.453  ym  band  and  the  0.745-0.785  ym  band  to  the 
existing  two  ocean  channels  of  the  scanning  radiometer  provides  the  quanti¬ 
tative  measurement  capability  of  chlorophyll  content  and  sediment  level  in 
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the  ocean.  Therefore,  the  improved  10-band  scanning  radiometer  will  have  an 
Integral  capability  of  monitoring  meteorological  phenomena,  water  color  as 
well  as  natural  disasters  such  as  floods,  fires,  and  drought. 

The  author  wishes  to  express  his  thanks  to  Mr  Kuang  Dingpo  and  Mr  Zheng  Qinpo 
for  their  assistance  in  this  work. 

Table  5.  Performance  of  the  10-Band  Scanning  Radiometer 


Channel  j 
No  1 

Band 

Spatial 

resolution 

Coverage  1 
period 

Application 

i 

1  i 

0.433-0.453  ym 

1.08x1.08  km2 

24  hours 

Ocean  chlorophyll 

2 

0.48-0.53  ym 

1.08x1.08  km2 

24  hours 

Ocean  chlorophyll 

3  ' 

0.53-0.58  ym 

1.08x1.08  km2 

24  hours 

Ocean  chlorophyll,  yellow 
substance 

4 

0.58-0.68  ym 

0.54x0.54  km^  j 

24  hours 

j  • 

Cloud  maps,  ground  cover, 
ocean  sediment 

5 

0.725-1.1  ym 

0.54x0.54  km2  ' 

;  24  hours 

Cloud  maps,  ground  cover, 
water  and  land  separation 

6 

0.745-0.785  ym 

1.08x1.08  km2 

j  24  hours 

Atmospheric  correction 

7 

1.58-1.64  ym 

1.08x1.08  km2 

j  24  hours  | 

Drought  detection,  cloud 
and  snow  separation 

8 

3.55-3.93  ym 

1.08x1.08  km2 

1  12  hours 

Fire  detection,  ocean 
surface  temperature 

9 

10.5-11.5  ym 

1.08x1.08  km2 

12  hours 

Infrared  cloud  maps, 
temperature  measurement 

10 

11.5-12.5  ym 

1.08x1.08  km2 

12  hours 

Infrared  cloud  maps, 
temperature  measurement 
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AEROSPACE 


Dynamics  of  Simplified  Model  of  Solar  Panel  and  Supporting  Components 

90FE0179B  Beijing  YUHANG  XUEBAO  [JOURNAL  OF  CHINESE  SOCIETY  OF  ASTRONAUTICS] 
in  Chinese  No  2,  30  Apr  90  pp  57-64 

[Article  by  Zhang  Jiazhong  [1728  0857  6988]  and  Shao  Chengxun  [6730  2052 
8133] ,  Harbin  Poly technical  University] 

[Text]  Abstract 

Experimental  results  and  theoretical  analysis  of  a  simplified  model  of  a 
solar  panel  and  its  supporting  components  used  on'  a  space  station  are 
presented.  Specifically,  the  dynamic  characteristics  of  the  system  when  it 
is  subject  to  disturbances  such  as  attitude  adjustment  and  orbit  maneuvers 
are  investigated.  The  results  show  that  because  of  the  inertial  non-linear 
coupling  effect  of  a  large  flexible  space  structure,  the  presence  of 
internal  resonance  will  produce  abrupt  changes  in  the  dynamic  response  of 
the  structure  when  the  amplitude  of  the  disturbance  torque  exceeds  a 
threshold  value.  These  changes  exhibit  non-linear  characteristics  of 
simultaneous  excitation  of  the  coupled  modes .  The  abrupt  increases  in  the 
amplitudes  are  expected  to  cause  difficulties  in  maintaining  the  shape  and 
controlling  vibration  of  the  flexible  structure. 

I .  Introduction 

For  a  large  space  structure  such  as  the  space  station  which  operates  in 
near-zero  load  conditions,  the  strength  requirements  are  very  low.  To 
reduce  the  weight  of  the  structure,  many  of  its  components  are  extremely 
flexible.  However,  these  components  are  subject  to  large  overload  during 
orbit  injection;  therefore  they  are  compactly  stowed  during  launch,  then 
deployed  or  re— assembled  after  achieving  orbit.  A  typical  example  of  such 
flexible  component  is  the  foldable  solar  panel  discussed  in  reference  1. 
Figure  1  shows  the  configuration  of  a  postulated  space  station,  which  con¬ 
sists  of  the  extremely  flexible  solar  arrays,  the  radiation  panels  and  a 
large  antenna. 

Because  a  space  station  generally  undergoes  a  series  of  on-orbit  operations 
such  as  deployment  and  re-assembly  of  its  components,  docking  with  the 
spacecraft,  attitude  control  and  orbit  maneuvers,  it  is  inevitably  subject 
to  different  types  of  disturbances.  A  large  flexible  structure  such  as  the 
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space  station  has  a  high  concentration  of  low-frequency  natural  modes.  An 
analysis  of  a  postulated  space  station^  shows  that  there  are  14  natural  modes 
in  the  frequency  range  0.21-0.75  Hz,  and  20  natural  modes  within  a  2-Hz 
range.  These  include  the  modes  of  the  flexible  structure  as  well  as  the 
aggregate  modes  of  the  space  station.  Furthermore,  these  flexible  components 
which  are  connected  through  the  main  body  of  the  space  station  are  mutually 
coupled.  The  non-linear  coupling  has  a  significant  effect  on  the  dynamic 
characteristics  of  the  system.  Reference  3  pointed  out  that  these  non-linear 
coupling  effects  will  cause  local  vibrations  in  the  structure.  References 
4  and  5  have  analyzed  the  effect  of  inertial  non-linear  coupling  on  the 
dynamic  characteristics  of  a  flexible  structure,  and  pointed  out  that  the 
presence  of  internal  resonance  in  the  system  will  cause  abrupt  changes  in 
the  system  dynamic  response  even  if  the  coupling  is  weak.  For  a  large 
flexible  structure  (Figure  1)  with  a  large  number  Of  natural  modes,  the 
presence  of  internal  resonance  is  almost  unavoidable.  Therefore,  this 
effect  must  be  carefully  studied. 


Figure  1.  Proposed  Space  Station 

The  dynamic  equations  of  a  large  space  structure  expressed  in  modal  coordi 
nates  are  of  the  following  form: 


=  +  -f2carSrgr 

+  e*[3rd  order  terms]  (1) 

where  e  is  a  small  parameter,  and  the  second  term  on  the  right  side  of  the 
equation  is  the  second-order  inertial  coupling  term.  Qj-  is  the  generalized 
force,  Mr  is  the  generalized  mass,  (Oj-  is  the  frequency  of  the  r-th  order 
linear  mode,  and  Arij  and  Brij  are  the  corresponding  coupling  coefficients. 
Under  normal  conditions,  the  effect  of  these  non-linear  coupling  terms  can 
be  neglected,  which  leads  to  the  conventional  linear  dynamic  equations: 


$r+2ar£r  jr  +  Orqr 


Qr 

"Mr 


(2) 
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Equation  (2)  provides  an  accurate  description  of  the  system  characteristics 
and  yields  the  correct  system  response^  But  when  internal  resonance  exists 
between  the  natural  modes  of  the  system,  i.e.,  when  the  relationship 

«(— O;  (3) 

exists  between  two  natural  modes,  or  when  the  relationship 

cok—cot+'an  (4) 

exists  between  three  natural  modes,  or  when  both  relationships  are  present, 
then  as  the  amplitude  of  the  disturbance  increases,  abrupt  changes  of  the 
system  dynamic  response  may  occur.  In  this  case,  the  effect  of  non-linear 
coupling  can  no  longer  be  neglected.  This  article  presents  the  experi¬ 
mental  results  and  theoretical  analysis  of  the  phenomenon  of  internal 
resonance  under  non-linear  coupling  for  a  simplified  model  of  solar  panel 
and  its  supporting  components. 

II.  Experiments 


In  order  to  study  the  effect  of  inertial  non-linear  coupling  on  system 
dynamics,  we  have  conducted  an  experimental  Investigation  and  performed 
theoretical  analysis  of  a  simple  model  of  solar  panel  and  its  supporting 
system  as  shown  in  Figure  2.  The  solar  panel  is  modeled  as  a  beam  under 
bending  and  torsional  loads,  the  supporting  system  is  modeled  as  another 
beam.  We  only  consider  the  case  where  the  solar  panel  is  perpendicular  to 
the  supporting  beam;  the  physical  parameters  of  the  two  components  are 
adjusted  so  that  the  natural  bending  mode  in  the  plane  of  the  solar  panel 
and  the  supporting  components  is  wi  =  6.1  Hz,  the  natural  torsional  mode 
about  the  axis  of  the  supporting  beam  is  wx  =  56  Hz,  and  the  natural 
torsional  mode  of  the  solar  panel  about  its  own  axis  is  =  52  Hz.  Clearly, 
these  three  frequencies  satisfy  the  following  internal  resonance  relation¬ 
ship  : 

Ot --Oi'f  COp  /c\ 
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A  sine-wave  excitation  signal  is  generated  using  a  SC0200  scan-controlled 
signal  generator  through  a  LDS  resonator;  its  spectrum  and  time  history  are 
recorded  using  a  spectrum  analyzer  and  a  X-Y  recorder.  The  torsional  vibra¬ 
tion  of  the  supporting  beam  is  measured  by  a  micro-accelerometer  located  at 
the  end  of  the  supporting  beam.  The  torsion  and  bending  of  the  solar  panel 
are  measured  by  the  strain  gauge  located  on  both  sides  of  the  solar  panel. 

Figure  3  shows  the  torsional  frequency  response  of  the  supporting  beam  under 
a  sine-wave  scan.  Figure  (a)  shows  the  response  under  low  excitation  level, 
which  is  clearly  linear  and  has  a  peak  response  at  56  Hz.  Figures  (b) ,  (c) 
show  that  as  the  excitation  level  increases,  the  non-linear  characteristics 
become  increasingly  apparent,  and  multiple  peaks  appear  on  the  frequency 
response  curve.  Figure  (d)  shows  the  typical  non-linear  behavior  of  the 
frequency  response  as  the  scan  direction  is  reversed. 


Figure  3 

Figure  4  shows  the  torsional  spectrum  of  the  supporting  beam,  and  the  bending 
spectrum  and  torsional  spectrum  of  the  solar  panel  under  a  56-Hz  sine-wave 
excitation.  It  can  be  seen  that  at  low  excitation  levels  (Figure  (a)),  only 
torsional  vibration  is  produced;  at  high  excitation  levels  (Figure  (b)),  a 
drive  frequency  of  w-p  will  excite  not  only  torsional  vibration  of  supporting 
beam,  but  also  coupled  bending  and  torsional  vibrations  of  the  solar  panel. 

In  other  words,  at  high  excitation  levels,  the  non-linear  effects  become  very 
pronounced;  the  torsional  disturbances  acting  on  the  supporting  beam  are 
transmitted  to  the  solar  panel  through  internal  resonance,  producing  bending 
and  torsional  vibrations  in  the  panel. 
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Figure  4 


Figure  5  illustrates  the  effect  of  the  disturbance  level;  curve  1  is  the 
torsional  vibration  of  the  supporting  beam,  and  curve  2  is  the  vibration  of 
the  solar  panel;  the  horizontal  axis  denotes  the  excitation  level.  These 
curves  clearly  show  that  when  the  excitation  level  is  less  than  a  certain 
critical  value,  the  system  dynamics  remain  in  the  linear  domain,  where 
torsional  vibration  is  produced  only  in  the  supporting  beam,  the  solar  panel 
is  free  of  vibration.  When  the  excitation  level  is  greater  than  this 
critical  value,  the  system  dynamics  become  non-linear;  the  disturbance  energy 
acting  on  the  supporting  beam  is  transmitted  to  the  solar  panel  through 
internal  resonance,  and  produces  vibration  in  the  panel.  As  the  excitation 
level  further  increases,  the  amplitude  of  the  torsional  vibration  of  the 
supporting  system  reaches  a  saturation  level,  whereas  the  amplitude  of  vibra¬ 
tion  of  the  solar  panel  continues  to  increase. 


Curve  1  denotes  the  torsional 


vibration  about  the  beam  1  axis; 


Curve  2  denotes  the  bending 
vibration  in  the  plane 


Figure  5 


For  a  large  flexible  structure  with  ntimerous  low-frequency  natural  modes, 
this  Internal  resonance  phenomenon  clearly  is  a  potentially  dangerous 
problem. 
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III.  Theoretical  Analysis 

The  solar  panel  and  its  supporting  system  are  represented  by  the  simplified 
model  shown  in  Figure  6.  Assume  that  the  disturbance  torque  is  about  the 
axis  of  the  supporting  beam  and  has  a  sine-wave  excitation  with  frequency  fJ. 
Let  the  supporting  system  and  the  solar  panel  be  represented  by  beam  1  and 
beam  2  respectively*  both  of  which  are  subject  to  bending  and  torsion. 


The  relative  deformation  of  beam  2,  denoted  by  u(z,t)  and  (p(z,t),  can  be 
expressed  as  follows: 

uiz,t) -Uoftis)  (6) 

(7) 


where 


Uo  =  h',  Kol 


(8) 


.  is  a  component  of  the  eigenvector  of  the  first-order  bending  mode  of  the 
system*  and  is  the  generalized  coordinate  of  the  bending  mode  in  the  plane. 
ijfQ  is  the  coordinate  of  the  torsional  mode  about  the  beam  2  axis.  Let  Oq  be 
the  generalized  coordinate  of  the  system  torsional  mode  about  the  beam  1  axis* 
then  the  system  dynamic  equations  expressed  in  modal  coordinates  are: 


«o— (afco— cosQt  —  2dra)rc()~  +  +  i<PoD 

J  0  0 


(9) 


if-'o  CO%^o —  ZSipCOv^po 


^  I  1  —  —  ZdiCOi^V 


-Ms. 

io 

M, 

■"  'M\ 


i  ipot 


(10) 

(11) 


where  is  the  second-order  coupling  coefficient: 

Bir  =  \‘  {l-z)f,{z)g{z)dz 
Jo 


(12) 
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It  can  be  seen  from  equations  (9) -(11)  that  the  coupling  term  + 

Wi^Po]  on  the  right  side  of  equation  (9)  absorbs  part  of  the  direct  excitation 
of  the  ao  response;  the  coupling  terms  and  oiofo  on  the  right  side  of 

equation  (10)  and  equation  (11)  cause  indirect  excitation  of  the  torsional 
deformation  about  the  beam  2  axis  nnd  of  the  bending  deformation  W-j_ 
respectively.  These  coupling  terms  are  important  if  the  following  internal 
resonance  condition  exists: 


CO  r  I  "t"  (13) 

Otherwise  they  can  be  neglected.  This  conclusion  and  the  experimental  results 
presented  in  the  previous  section  can  be  verified  by  theoretical  analysis  of 
this  simplified  model  using  the  multiple  scale  method. 

To  proceed,  we  introduce  an  arbitrary  small  parameter  e  and  re-write 
equations  (9)-(ll)  in  the  form: 


Y  =  sLrcosQt-2iTC0Tf  1-  l^kniX  iX  ^-^2X  ,  +  X  ,X  2)1 

X,  +  a)f;r.=£[-2|.o,X,-A'nF:.Y;] 

2  "f"  coV  ~  £  C  ^^12^^  i-Yi!l 

where 

eY,—ao>  £/Ai=W  1,  ElX:  =  rg:s, 

M  —  I  m  A;:  =  I-Sic 

=  £5  =  <5 

Jo  ''0 


(14) 

(15) 

(16) 


(17) 


r  is  the  radius  of  gyration,  is  the  generalized  mass  corresponding  to  the 
system  bending  mode,  Ig  is  the  polar  moment  of  inertia  of  beam  2,  Mg  is  the 
mass  at  the  end  of  beam  2,  Tg  is  the  amplitude  of  the  applied  torque,  Jg  is 
the  generalized  rotational  inertia,  yo  is  a  component  of  the  torsional  mode 
at  the  point  of  the  applied  torque,  and  I  is  the  length  of  beam  2. 

By  Introducing  a  slowly  time-varying  parameter  Ti  =  et,  the  solution  of 
equations  (14)-(16)  can  be  expressed  in  the  form: 

F,  it)  =6.  (T,)  cosLcort  +  ^:  (T,)  ]  +eF.,  (r,r.) 

-Y.  (0  =  fli  (r.)  cos  [o.f  +  a,  (r, )  ]  -r  e  A', ,  (t,  T  ■, ) 

Xi{i)  =02(7 ,)  cosLcpvt  +  OjiT ,)]  +  eA’',,(f, 7,) 

Equations  (18) -(20)  represent  the  torsion  of  beam  1,  the  bending  of  the 
system  in  the  plane,  and  the  torsion  of  beam  2  respectively.  It  can  be  seen 
from  these  equations  that  each  solution  corresponds  to  a  domain  of  its  own 
natural  frequency.  From  the  bounded  conditions  of  the  disturbance  term 
Yi]^(t,  Tj^) ,  one  can  derive  a  first-order  differential  equation  for  each  of 
the  response  amplitudes  b^,  ai  and  a2  with  respect  to  the  slowly  varying  time 
T^.  The  degree  of  resonance  can  be  defined  by  two  parameters  a  and  02’- 


(18) 

(19) 

(20) 


Q  =  cot^£0’ 

(21) 

C0T  —  C0t’^(fv  +  £cri 

(22) 

Expanding  the 
one  obtains: 

differential  equations  for  the  modal  responses  ai,  ^2 

,  and  hi , 

coia'i  =  -  ^tQ)]ai-\-  ^^Ki:j)Ta2bzSiny2 

(23) 

Qyiata[  —  iioyraibzOOsVz 

(24) 

o<pai~  ^<pco\>ai'^ K a ib 2  siny2 

(25) 

Otpaicct  “  j cosy^ 

(26) 

0) r6l  ^  /  6  — /C , 2  (cy  1  -  0)^ )  "ai 02  sin  y ,  -h  sm  y j 

(27) 

coTbifi\  — — ^  u K 12  {co (^^<p)  “ Gicizcosy 2  2 

(28) 

where 

Yi-aTi-jS, 

V2  =  cj,T,— 

(29) 

(30) 

In  the  above  equation,  the  superscript  denotes  derivative  with 

The  system  equations  (23) -(28)  have  two  steady-state  solutions 
the  solutions  is: 

respect  to 
.  One  of 

ai  =  0,  a!  =  0 

(31) 

,  ^ . 

(32) 

which  indicates  that  the  indirectly  excited  motions  are  stationary,  and  the 
response  of  the  directly  excited  motions  is  linear. 

The  second  steady-state  solution  is: 


(33) 

(34) 

(35) 
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where 


jCl!®r(liOi  +  S ?£<)»’)  ~o{g  +  (7,)] 
+|fO,.)CT  +  |rOr(o-  +  cr,)3 

y--  ^ _ ScOiOa-y/ Slip 

fj  12  (s  1®  1 


(36) 

(37) 

(38) 


This  solution  represents  the  non-linear  response  of  the  system.  The  three 
types  of  motions  are  excited  simultaneously;  also,  when  a  =  a^,  the 
torsional  amplitude  has  a  minimum  value  which  is  independent  of  the  level 
of  excitation.  This  is  the  so-called  saturation  phenomenon  which  is  clearly 
evident  in  the  experimental  results  of  Figure  5.  The  amplitudes  of  a^  and  a2 
increase  with  increasing  excitation  level.  This  phenomenon  can  also  be 
observed  in  Figure  5. 

Equation  (35)  also  shows  that  whether  or  not  a  stable  solution  exists  depends 
on  the  system  parameters. 


The  To  value  of  the  critical  point  can  be  obtained  from  equations  (32)  and 
(33) : 


To  = 


8coi<j}<p  j 
COt  y 


(a^  +  iWr) 


(39) 


Clearly,  To  is  a  minimum  when  +  a  =  0  and  a  «  0: 


T'om!- 


(40) 


However,  +  a  =  0  and  a  =  0  imply  that  the  Internal  resonance  relationship 


iV  r  —  O);  “f  co<p 


and  the  resonance  relationship  0  =  0).^  exist. 

Figure  7  shows  the  theoretical  curves  of  the  amplitudes  of  the  modal 
responses  as  a  function  of  the  excitation  parameter  t.  When  t  exceeds  the 
critical  value,  the  amplitudes  of  the  coupled  modal  responses  vary 
differently:  one  approaches  saturation,  and  the  other  two  undergo  abrupt 
changes,  then  continue  to  rise  sharply. 

IV.  Conclusion 

In  this  article,  the  effect  of  Inertial  non-linear  coupling  on  the  dynamic 
characteristics  of  a  flexible  system  is  studied  using  a  simplified  model  of 
solar  panel  and  its  supporting  system.  Experimental  results  and  theoretical 
analysis  show  that  when  the  external  excitation  is  in  resonance  with  one  of 
the  non-linear  coupled  modes,  and  internal  resonance  also  exists  between  the 
coupled  modes,  the  non-linear  coupling  can  produce  abrupt  changes  in  the 
d3niamic  characteristics  of  the  system.  The  effect  of  internal  resonance 
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produces  indirect  excitation  of  the  coupled  modes  which  are  independent  of 
the  external  excitation.  For  a  large,  flexible  structure  such  as  the  space 
station,  there  are  numerous  sources  of  excitation  and  a  high  concentration  of 
natural  modes  which  may  cause  internal  resonance.  This  type  of  internal 
resonance  is  a  potential  threat  to  vibration  suppression,  configuration 
maintenance,  pointing  accuracy  and  control  accuracy  for  a  large  flexible 
structure.  Therefore,  detailed  investigation  of  this  effect  of  non-linear 
coupled  vibration  is  essential  to  the  design  of  space  station. 


Excitation  parameter  t 

Figure  7 
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[Text]  Abstract 

A  method  of  coarse  geometric  rectification  of  full-size  satellite  panoramic 
images  is  proposed.  The  method  is  based  on  certain  intrinsic  geometric 
relationships  between  the  satellite  image  and  the  ground,  and  takes  into 
account  the  physical  process  of  satellite  imaging.  The  rectification  of 
satellite  images  is  accomplished  by  applying  parameter  optimization  technique 
which  does  not  depend  solely  on  the  satellite  parameters,  but  also  relies 
on  a  small  number  of  carefully  selected  and  accurately  measured  ground 
control  points  as  well  as  a  geometric  model.  Simulation  results  show  that 
the  proposed  method  is  very  effective  for  the  coarse  processing  of  satellite 
panoramic  images. 

I .  Introduction 

Ground  observation  and  imaging  by  a  retrievable  remote-sensing  satellite  is 
accomplished  by  panoramic  scanning.  Because  of  the  low  satellite  orbit  and 
the  large  scan  angle,  the  images  suffer  severe  geometric  distortions, 
particularly  at  the  edges  where  errors  can  be  as  large  as  several 
kilometers.  In  order  to  promote  and  apply  the  satellite  panoramic  imaging 
technique,  it  is  essential  to  find  a  practical  method  of  geometric  rectifi¬ 
cation. 

The  geometric  methods  which  are  currently  used  for  satellite  image 
processing  include  the  model  method  and  the  polynomial  method. ^  The  model 
method  involves  establishing  a  mathematical  model  to  describe  the  satellite 
imaging  process  and  requires  accurate  knowledge  of  the  satellite  orbit, 
satellite  attitude  and  camera  parameters;  the  accuracy  of  these  parameters 
will  have  a  significant  effect  on  the  accuracy  of  geometric  processing.  The 
polynomial  method  does  not  require  knowledge  of  satellite  parameters,  but  it 
does  require  a  large  number  of  accurately  measured  and  uniformly  distributed 
ground  control  points;  this  method  produces  rather  poor  results  at  the  edges 
of  the  images. 
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The  method  proposed  in  this  paper  is  based  on  the  geometric  relationships  of 
satellite  imaging  and  tries  to  avoid  the  deficiencies  of  the  above  two 
methods.  Specifically,  it  takes  into  account  the  physical  process  of 
satellite  imaging,  but  does  not  rely  solely  on  the  satellite  parameters; 
rather,  it  uses  the  intrinsic  geometric  relationship  between  the  satellite 
image  and  a  small  number  of  accurately  measured  ground  control  points  and 
performs  geometric  rectification  of  the  satellite  image  by  applying  optimal 
parameter  estimation  technique. 

II.  Mathematical  Model^ 

To  simplify  the  analysis,  the  earth  is  assumed  to  be  a  sphere  of  radius  R, 
and  the  altitude  of  the  satellite  is  denoted  by  H.  Since  the  effective  area 
of  each  film  is  approximately  200x1400  mm2,  its  width  corresponds  to  a  size 
of  only  several  tens  of  kilometers  on  the  ground.  Hence,  the  projected  scan 
lines  can  be  approximated  by  straight  lines,  and  the  ground  segment  corres¬ 
ponding  to  each  film  can  be  approximated  by  a  curved  surface  taken  from  a 
cylinder  of  radius  R,  as  indicated  by  the  surface  EFPB  in  Figure  1.  Clearly, 
for  a  scan  angle  S,  one  can  apply  the  law  of  sine  to  the  triangle  AOSPq  to 
obtain  the  following  equation: 

2a=  sin‘'|^sin  • /l+  (1) 

■  R  ■  sin2g 

sin^  ' 

Recognizing  that  AOSiPo  is  an  isosceles  triangle,  one  can  determine  the 
corresponding  angles,  and  apply  the  law  of  sine  to  the  triangle  ASSj^Pq  to 
obtain: 

1  =  H  cosa/cos(a  + P)  (3) 


Figure  1.  Geometric  Relationship  of  Satellite  Imaging 
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1,  Geometric  Relationship  Along  the  Scan  Direction 

As  shown  in  Figure  2,  in  the  ground  coordinate  system  0-XY,  Lq  is  the  ground 
projection  of  the  panoramic  scan  slit  when  the  scan  angle  g  =  0;  it  is  repre¬ 
sented  by  the  equation  Y  =  KX+E.  At  time  t^,  the  ground  projection  of  the 
slit  is  L±»  and  an  image  of  the  straight  line  L^*  which  passes  through  the 
point  (Xi,  Y^)  is  taken.  The  distance  between  and  Lq  is  vt^,  where  v  is 
the  translational  velocity  of  Lq  due  to  earth  rotation,  which  is  to  be  deter¬ 
mined.  It  is  a  function  of  the  inclination  of  the  satellite  orbits  the 
satellite  velocity,  the  earth  rotational  speed  and  the  latitude  and  longitude 
of  the  satellite  subpoint.  Dj[  is  the  distance  between  Lf’  and  Lj[. 


r 


Figure  2.  Geometric  Relationship  in  the  Ground  Coordinate  System 

In  the  original  film  coordinate  system  0-xy  shown  in  Figure  3,  it  is  assumed 
that  the  satellite  attitude  remains  unchanged  during  the  imaging  process;  in 
other  words,  during  this  period  the  equation  of  the  theoretical  center  edge- 
contact  line  (3=0)  remains  constant  Aq:  y  =  kx+e.  By  "theoretical"  it  is 
meant  that  Hq  passes  through  the  primary  image  point,  and  the  image  line  at 
any  time  is  parallel  to  Aq.  df  is  the  distance  between  and  Aq,  and 
corresponds  to  the  ground  distance  D^. 


Figure  3.  Geometric  Relationship  in  the  Film  Coordinate  System 
From  the  above  analysis,  one  can  readily  derive: 


\KX,  +  E-Y<\ 
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c/,-= 


\kxi  +  e  —  yi\ 


■v/  1  + 


(5) 


where  K,  E,  k,  e,  v  are  to  be  determined,  t^  »  (Vq  /  1+k^) ,  Vq  is  the  film 
Speed  of  the  camera. 

Substituting  =  di/f  (f  is  the  focal  length  of  the  camera)  into  equation 
(1),  one  can  derive  the  following  relationship  between  and  d^: 


£),  =  i?  .  {  sin-[sin-y''(l+  y)]- 


(6) 


It  is  clear  that  even  if  errors  are  present  in  the  satellite  attitude  and 
orbit  parameters,  as  long  as  they  remain  unchanged  during  the  imaging 
process,  this  geometric  relationship  also  remains  unchanged.  In  other  words, 
equation  (6)  shows  that  the  accuracy  of  geometric  processing  is  independent 
of  errors  in  the  satellite  attitude  parameters,  which  is  a  major  advantage 
of  this  method.  Since  each  film  only  takes  3  seconds  to  complete  during  the 
imaging  process,  it  is  reasonable  to  assume  that  the  effective  parameters 
remain  unchanged  during  this  period. 

2.  Geometric  Relationship  Perpendicular  to  the  Scan  Direction 

To  account  for  small  variations  in  satellite  motion,  earth  rotation  and  the 
direction  of  film  advancement,  we  use  a  second-order  curve  to  represent  the 
scan  line  of  the  slit-center  and  its  ground  track;  these  equations  are  of 
the  form:  x  =  ay^fby+c  (film)  and  Y  =  AX2+BX+C  (ground) ,  as  Indicated  in 
Figure  2  and  Figure  3.  In  the  ground  coordinate  system,  let  us  draw  a 
straight  line  from  the  ground  control  point  (Xi,  Yi):  (Y-Yi)  =  K(X-Xi) , 
which  intersects  the  center  scan  line  Y  “  AX^BX+C  at  (Xq*  Yq)  ,  then  we  have: 

DYr^^/'CX,^oyTiV7-Y'y  (7) 

Similarly,  in  the  original  film  coordinate  system,  we  draw  a  straight  line 
from  (xi,  yi) :  y-yi  =  k(x-xi) ,  which  Intersects  the  center  scan  line 
X  =  ay2+by+c  at  (Xq,  yo) .  then  we  have: 


dyt-V  (xi-xo)^~  {yi-y»)^ 


(8) 


Now,  we  shall  derive  a  geometric  relationship  between  DYf  and  dy^.  Figure  4 
shows  the  geometry  between  DYi  and  dy^  where  a  pitch  angle  r  is  present. 

0  is  the  angle  of  the  scan  beam.  It  should  be  pointed  out  that  since  ^  and 
dy  <-aT>  be  either  positive  or  negative,  they  can  be  combined  to  form  several 
different  geometries;  however,  the  derivation  in  this  paper  is  based  only  on 
the  geometry  shown  in  Figure  4;  it  can  be  readily  extended  to  other 
geometries.  From  Figure  4,  it  is  clear  that 
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(9) 


DY  =  l  ■  [tg(95  +  0)  -  tg^)] 

=  p .  r + 
sin/9  Li— tg9;tg0 

=  M  .  _ '^0?" . T  f'^S<p-rdy  ,  -f 

cos  a+j-^f-^S<fdy 

Equation  (9)  is  affected  only  by  the  pitch  angle  4';  it  is  independent  of  other 
attitude  parameters. 


Figure  4.  Geometric  Relationship  Between  DY  and  dy 
3.  Objective  Function 

In  this  paper,  the  physical  parameters  are  determined  by  solving  the 
parameter  optimization  problem  using  the  Powell  method.  The  objective  func¬ 
tions  to  be  optimized  are: 


■  =  '  I 


cos  a 


4.  Geometric  Rectification 


r.  fJS£y:dyj 

•  dy. 


(10) 


(11) 


Having  determined  the  parameters  {K,  E,  k,  e,  v.  A,  B,  C,  a,  b,  c,  «»’}  using 
the  Powell  method,  the  process  can  be  reversed  to  establish  the  relation¬ 
ships  between  corresponding  points  on  the  original  film  and  the  new  film 
after  rectification.  A  difficult  problem  with  this  process  is  establishing 
the  quadrant  relationship  between  corresponding  points  in  the  new  picture  and 
the  original  picture. 
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Consider  a  point  p(x,y)  on  the  original  film  and  its  corresponding  point  on 
the  new  film  P(X,Y).  The  computational  procedure  and  the  key  formulas  are 
given  below: 

1)  Determine  Xq,  yo 


j(>'-Vo)/  r -.Vo) 
1  Vo  =  -  o  \'o  +  c 


1  >  I’o) 


2)  Calculate  d,  dy 


i _ \kx^  e  —  y\ 

1  +  A' 

dy  =  V'x-x^)-^  {y ->’■>)- 


3)  Calculate  D,  DY  from  equations  (6)  and  (9) 

4)  Solve  the  simultaneous  equations 

Y-Yo^K{X-X,) 

Y,^AXli-BK,  +  C 

DY  =V{X-X,)'--r{Y-Y,)-  ^  ( A\, .  I'o ,  -Y ,  Y) 

\KX  +  E-Y\__ 

This  completes  the  transformation  between  the  points  p(x,y)  and  P(X,Y).  It 
should  be  noted  that  the  positional  relationship  between  corresponding  points 
on  the  original  film  and  the  new  film  is  controlled  by  the  four  "quadrants" 
formed  by  the  center  edge  contact  and  the  center  scan  line;  i.e.,  the 
"quadrant"  relationship  between  corresponding  points  remains  unchanged. 

To  save  computer  time  in  the  actual  rectification  process,  we  only  use  the 
above  model  to  calculate  the  coarse  lattice  points;  then  we  apply  the  finite 
element  method  (first-order  approximation  inside  a  triangle)  to  calculate 
points  within  a  lattice. 

5.  Other  Considerations 

Since  the  panoramic  scan  camera  uses  mechanical  drive  mechanisms  for  film 
advancement,  jitters  will  inevitably  appear,  causing  errors  in  the  images. 

To  solve  this  problem,  we  introduce  in  the  objective  function  an  additional 
parameter  6,  which  is  a  function  of  d^.  The  form  of  this  function  is  shown 
in  Figure  5,  where  A  is  the  residue  at  d^,  and  L  is  the  range  of  influence 
of  jitter,  which  is  the  parameter  to  be  optimized.  For  a  control  point  whose 
residue  is  larger  than  a  certain  value,  optimization  is  performed  on  an 
objective  function  with  the  added  parameter  6(di),  which  determines  a  model 
function  that  reflects  the  jitter  of  the  entire  frame. 
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Figure  5.  Simulated  Jitter  Function  of  Film  Advancement 

In  addition,  because  of  variations  in  satellite  attitude  parameters  during 
the  scan  process  and  because  of  dimensional  considerations,  the  film 
digitization  process  must  be  carried  out  frame  by  frame  (due  to  the  limita¬ 
tion  of  C-4500) .  Therefore,  the  center  scan  line  or  its  ground  track  cannot 
be  modeled  by  a  single  second-order  curve;  it  must  be  modeled  within  each 
frame  to  achieve  the  required  processing  accuracy. 

Another  factor  that  should  be  considered  is  the  effect  of  difference  in 
surface  altitude.  This  effect  is  particularly  pronounced  for  domestic 
satellites  as  they  are  in  low-earth  orbits  and  have  large  scan  angles.  In 
this  article,  the  effect  of  difference  in  surface  altitude  is  taken  into 
account  by  introducing  a  correction  factor  to  the  orbit  altitude  H  in  the 
geometric  relationship. 

III.  Test  Results 

In  order  to  demonstrate  the  feasibility  of  the  proposed  method,  we  have 
conducted  test  on  a  full-size  satellite  panoramic  image.  The  dimensions  of 
the  image  are  approximately  200x1400  mm2.  The  image  is  divided  into  six 
frames,  each  of  which  is  digitized  using  the  C-4500;  the  individual  frames 
are  joined  together  in  the  computer  using  Identical  points  to  form  a  com¬ 
plete  digitized  image.  The  control  points  are  selected  by  performing  a 
correlation  analysis  between  the  image  and  the  terrain  map.  However, 
because  a  large  portion  of  the  area  covered  by  this  image  is  desert, 
selection  of  control  points  is  very  difficult,  and  they  are  not  uniformly 
distributed.  As  a  result,  the  RMS  processing  errors  using  the  polynomial 
method  or  the  traditional  model  method  are  of  the  order  of  800-1000  m, 
which  clearly  do  not  meet  the  accuracy  requirements.  In  this  paper,  a  total 
of  50  ground  target  points  are  selected,  30  of  which  are  control  points,  and 
the  remaining  20  are  check  points.  After  processing  the  data  using  the 
method  described  in  this  paper,  the  RMS  error  is  a  =  206  m,  which  fully 
meets  the  accuracy  requirement  of  coarse  geometric  processing  for  satellite 
images . 

This  method  can  also  be  used  for  precision  rectification  of  sub-images  of  a 
region  with  a  large  number  of  uniformly  distributed  control  points.  For 
example,  a  test  was  performed  on  a  110x110  mm2  sub-region  using  30  control 
points;  the  RMS  error  of  the  rectified  image  is  a  =  45  m. 
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Higher  processing  accuracy  can  be  achieved  for  a  sub-image  because:  1)  there 
are  larger  number  of  control  points;  2)  the  various  parameters  remain  essen¬ 
tially  constant  within  a  shorter  time  interval;  and  3)  it  is  easier  to  model 
the  center  scan  line  and  the  center  edge  contact.  However,  precision 
processing  the  full-scale  image  all  at  once  would  be  difficult. 

IV.  Conclusion 

The  above  test  results  show  that  the  proposed  method  is  effective  for  coarse 
geometric  processing  of  satellite  panoramic  images.  This  method  is  particu¬ 
larly  useful  in  cases  where  the  model  method  is  not  applicable  because 
accurate  satellite  parameters  or  orbit  parameters  are  not  available,  or  where 
the  polynomial  method  is  not  applicable  because  uniformly  distributed  control 
points  are  difficult  to  obtain. 

Another  advantage  of  the  proposed  method  is  that  it  can  perform  simultaneous 
rectification  on  a  full-scale  200x1400  mm2  image  without  having  to  process 
the  individual  sub-images.  This  not  only  improves  processing  efficiency  but 
also  eliminates  possible  errors  introduced  by  the  joining  process  of  recti¬ 
fied  sub-images. 

It  is  generally  not  difficult  to  select  more  than  30  ground  target  points  on 
a  full-scale  image.  A  sufficiently  large  nvimber  of  uniformly  distributed 
ground  target  points  should  be  chosen  at  the  two  sides  of  a  satellite 
panoramic  image  where  large  distortions  occur.  Simulation  results  show  that 
this  method  can  meet  the  accuracy  requirements  for  geometric  processing 
satellite  images. 

The  authors  wish  to  express  thanks  to  Comrades  Ju  Tao,  Guo  Zhaozeng,  and  Li 
Lingzhi  for  their  support  in  developing  this  method. 
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[Text]  Abstract 

In  this  article,  an  Improved  method  of  extracting  the  poles  of  underground 
targets  is  presented.  It  is  pointed  out  that  in  a  noise  environment,  the 
eigenvalue  method  is  more  accurate  in  determining  the  number  of  poles;  it  is 
also  shown  that  the  eigenvalue  method  always  yields  smaller  error  of  approxi¬ 
mation  than  classical  methods.  Experiments  of  extracting  the  poles  of  buried 
targets  have  been  performed,  from  which  several  significant  conclusions  can 
be  drawn. 

I .  Introduction 

I 

An  effective  approach  for  detecting  underground  targets  is  the  transient 
electro-magnetic  field  method.  This  method  has  applications  in  many 
different  areas  such  as  hydrology,  geology,  resource  exploration,  detection 
of  buried  articles,  mine  detection,  and  anti-stealth  technology.  However, 
to  use  this  method,  it  is  necessary  to  solve  the  problem  of  target  identifi¬ 
cation  using  the  transient  target  response.  Since  the  poles  are  invariant 
properties  of  a  target,  they  can  be  used  as  effective  signatures  for  target 
identification. 1*2.3  in  1975,  Rai  Mittra  et  al.  proposed  the  Prony  method  of 
extracting  poles  from  transient  pulse  response.  This  method  attracted  a 
great  deal  of  interest,  but  further  study  showed  that  in  the  presence  of 
noise  and  measurement  errors,  the  accuracy  of  this  method  is  so  poor  that  it 
may  be  of  little  practical  use.  Therefore,  efforts  are  being  made  to  explore 
new  methods  of  pole  extraction. 

II.  Comparison  of  Square  Approximation  Error 

In  1978,  Van  Blaricum  et  al.  proposed  a  new  method  of  pole  extraction^  where 
the  problem  is  formulated  as  one  of  finding  the  smallest  eigenvalue  (Xmin) 
and  the  corresponding  eigenvector  a  of  the  matrix  H^H: 

aM//)  "(//)«  =  A  (1) 


26 


where 


satisfies  the  condition; 
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where  N  is  the  number  of  poles.  Once  is  determined,  the  poles  in  the 
Z-plane  can  be  obtained  from  equation  (2). 

Let  e  =  (to*  ei,  Sm)^  be  the  error  vector,  which  can  be  expressed  by  the 

following  relation: 


iH)a  =  £ 

If  W(e)  is  the  square  approximation  error,  then  we  have: 


(3) 


■  ,  •  ;  =  0 


(4) 


Different  a  vectors  may  yield  different  values  of  square  approximation  error; 
the  (X  vector  determined  by  the  eigenvalue  method  is  the  eigenvector  that 
corresponds  to  the  smallest  eigenvalue  Xm-in  of  the  matrix  rTh,  i.e., 

«’•(//) ’■(H)  a  =  Xr.ln  (5) 

Therefore,  the  problem  that  needs  to  be  investigated  is  the  relative  magni¬ 
tude  between  Xniin  W(e). 

By  arranging  the  eigenvalues  of  H^H  in  ascending  order: 

=  (6) 


then  it  is  always  possible  to  find  a  series  of  (N+1)  orthogonal  unit 

vectors  Ei  (i  =  0,  1,  - ,  N)  such  that  (u)  =  (Eq,  E^,  **•,  Eu)  is  an 

orthogonal  matrix,  and  satisfies  the  condition: 

0\ 

(m)’'(H)’'(H)(«)=  ♦... 

\  0  '  u 


(7) 


Consider  a  normalized  vector  a  =  (aQ,  a^,  •••,  from  matrix  theory  we 

can  write: 
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(8) 


H 

a  =  ^aiEi  or  a={u)a 

i  =-0 


Therefore, 


jris)=a^iHyiH)a 

=  ayu)ayH)yH)  (u)a 
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(9) 


It  can  be  seen  from  equation  (9)  that  when  all  the  eigenvalues  are  equal  to 
^min»  have  the  following  inequality: 

,v  m 

W  it)  (10) 

(  =  0  t  =  0 


Since  a  is  a  normalized  vector,  it  follows  that 

]rit)^^mi«  (11) 

Thus,  we  have  proved  that  the  square  approximation  error  of  the  eigenvalue 
method  is  always  smaller  than  that  of  classical  methods. 

III.  Determining  the  Number  of  Poles 

The  number  of  poles  can  be  easily  determined  using  the  orthogonality 
property  of  the  eigenvectors. 

1.  In  the  noise-free  case,  let 


then  the  following  procedure  can  be  used  to  determine  the  number  of  poles: 

a.  Construct  a  y  x  y  matrix  (<{>)  from  the  sample  values. 

b.  Find  the  eigenvalues  of  the  matrix  (i|))  and  determine  the  number  of  zero 
eigenvalues;  if  there  are  L  zero  eigenvalues,  then  the  number  of  poles  is 

V  =  y  -  L;  if  there  is  no  zero  eigenvalue,  then  the  dimension  of  (<J))  must  be 
increased  to  form  a  new  matrix  ((fi) . 

c.  When  there  is  only  one  zero  eigenvalue,  i.e.,  when  L  =  1,  then  the 
corresponding  eigenvector  is  the  Prony  vector  a,  and  the  poles  are  deter¬ 
mined  along  with  V. 

2.  In  the  case  where  noise  is  present,  this  method  can  still  be  applied. 

Let  the  noise  sequence  be  zero  mean  with  variance  since  the  noise 
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sequence  is  uncorrelated  in  time  and  also  uncorrelated  with  the  signal 
sequence,  the  orthogonality  property  of  the  eigenvectors  is  preserved. 

Analysis  has  shown  that  the  noise-free  matrix  (cji)  and  the  noise-present 
matrix  (({)')  are  related  by  the  following  expression: 

{(/>')  =  (<f>)  +Vcf^l  (12) 

where  I  is  the  unit  matrix.  Therefore,  the  noise-free  procedure  can  still  be 
used,  except  that  in  this  case,  the  eigenvalue  sought  is  not  zero.  Once 

an  eigenvalue  of  yo"^  is  found,  the  corresponding  eigenvector  is  the  a  vector. 

The  above  analysis  shows  that  this  method  is  very  effective  in  treating  the 
noise  problem;  once  the  two  basic  statistical  parameters  (mean  and  variance) 
of  the  noise  sequence  are  known,  the  number  of  poles  N  and  the  a  vector  can  be 
readily  determined. 

IV.  Experimental  Procedure  and  Data  Processing 
1.  Experimental  Set-Up  and  Procedure 

A  block  diagram  of  the  experimental  set-up  is  shown  in  Figure  1.  The  pulse 
width  is  T  =  1.2  ns,  the  amplitude  is  40  v,  and  the  repetition  frequency  is 
10  KHz;  the  number  of  sample  points  is  256,  the  time  duration  is  20  ns,  and 
the  sample  interval  is  0.078  ns.  The  environment  is  a  100  cm  x  100  cm  x  60  cm 
sand  box.  The  transient  pulses  are  transmitted  by  a  narrow-pulse  generator; 
the  transmitting  and  receiving  antennas  are  quasi-TEM  horn  antennas.  The  two 
receiving  antennas  have  identical  characteristics ;  the  received  signal  is 
amplified  and  stored  in  a  micro-computer.  The  synchronization  of  the  trans¬ 
mitted  pulse  and  the  sampling  circuit  is  also  controlled  by  the  micro¬ 
computer  . 


mitting  antenna 


antenna  antenna 

Figure  1.  Detection  of  Underground  Targets  Using  Transient  Pulses 
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2.  Processing  of  the  Initial  Waveform 

a.  Elimination  of  the  effect  of  ground  reflection.  In  this  experiment,  a 
dual  receiving  antenna  configuration  is  used  so  that  the  effect  of  ground 
reflection  can  be  eliminated  by  differencing  the  two  received  signals. 

b.  Noise  processing.  Noise  reduction  is  accomplished  by  averaging  several 
waveforms  measured  at  the  same  antenna  position.  Any  d.c.  component  is 
eliminated  by  adjusting  the  baseline. 

c.  Gate  switching.  Gate  switching  technique  is  used  to  eliminate  the 
direct  coupling  effects  of  the  antenna  and  premature  target  response. 

d.  Filtering  of  high-frequency  components.  Because  of  inconsistencies  in 
sampling  times  between  the  two  channels  of  the  received  signal,  high- 
frequency  "burrs"  often  appear.  To  solve  this  problem,  a  low-pass  filter 
is  used  to  process  the  waveform. 

The  above  pre-processing  of  the  data  significantly  improves  the  signal  to 
noise  ratio,  and  consequently  ensures  the  success  of  pole  extraction. 

V.  Conclusion 


Experiments  of  extracting  the  poles  of  various  targets  have  been  conducted; 
four  of  these  targets  are  shown  in  Figure  2.  The  targets  are  buried  in 
homogeneous  dry  sand  at  a  depth  of  3—10  cm;  the  receiving  and  transmitting 
antennas  are  separated  by  a  distance  of  40  cm,  and  the  sampling  interval  is 
0.078  ns.  The  measured  and  recovered  waveforms  are  shown  in  Figure  3,  where 
the  pole  locations  are  also  indicated. 


Target  1 
(plastic 
sphere) 


Target  3 
(plastic 
cylinder) 


Target  2 

(metal 

cylinder) 


Figure  2.  Physical  Properties  of  Four  Different  Targets 
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On  the  basis  of  analysis  and  experimental  results,  the  following  conclusions 
can  be  drawn: 

1.  The  eigenvalue  method  always  yields  smaller  error  of  approximation  than 
classical  methods. 

2.  In  the  presence  of  noise,  the  eigenvalue  method  can  still  be  used  to 
determine  the  number  of  poles  and  to  extract  the  poles  as  long  as  two  basic 
parameters  of  the  noise  statistics  are  known. 

3.  Changes  in  target  depths  (within  10  cm)  has  little  effect  on  the  poles j 
however,  the  measured  results  with  respect  to  target  depth  are  not  presented 
here  due  to  space  limitation. 

4.  The  initial  waveform  must  be  pre-processed  to  facilitate  pole  extraction. 

5.  The  initial  phases  of  conducting  and  non-conducting  targets  show  clearly 
opposite  sense  (see  Figure  3-a  and  Figure  3-b) . 

6.  It  can  be  seen  from  Figure  3  that  the  larger  the  scattering  object,  the 
smaller  the  real  part  of  the  pole  and  the  larger  the  imaginary  part}  con¬ 
versely,  the  smaller  the  scattering  object,  the  larger  the  real  part  and  the 
smaller  the  imaginary  part. 

The  authors  wish  to  thank  Mr  Chen  Jinxiong  and  Mr  Huang  Peikang  of  the  No  2 
Institute  of  the  Ministry  of  Aeronautics  and  Astronautics  Industry  for  their 
guidance  and  support;  they  also  wish  to  thank  Comrades  Cao  Chang jun,  Zi 
Haiyen,  Liang  Tao,  Li  Xun,  and  Fang  Y1  of  this  institute  for  their  assistance 
in  conducting  the  experiment. 
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Abstract  In  this  paper ,  a  pure  digital  control  method  for  the  outer  axis  of  a  three— axis 
flight  simulation  table  is  .worked  out  and  a  state  observer  is  used  to  obtain  the 
acceleration  signal  .The  method  is  simple  and  easy  to  be  emplayed,  and  especially  suitable 
for  the  electrohydrauUc  aervo— system  with  last  response  .  The  simulation  results  show 
that  the  system  given  in  this  paper  has  optimal  three— order  ^tem  performance  . 

Key  Wirdsc  Electrohydraulic  position  servo  system;  state  feedback;  digital  control . 

0  Intrcxiuction 

A  flight  simulation  table  is  a  kind  of  experimental  equipmoit  used  for  simulating  the  flying 
attityHff  of  atactical  missfle.  Therefoie^  it  is  required  to  have  the  petfonnanoe  of  high  accuracy, 
fast  response,  less  following  tiiw  and  little  overshoot  Now  a  days,  the  optimal  linearization 

tested  by  an  analog  potentiometer  can  only  reach  ±0.02%  (and  the  best  position  accuracy 

is  ±4*’32  O'Besides,  the  of  parameters  shifting  and  nonlinearity  often  occur  in 

odier  analog  oon^nents.  Therefore;  the  analog  oontibl  system  today  is  dif&cult  to  satisfy  the 
performance  required.  Furthermore;  because  the  dasskal  digital  control  is  only  restricted  to  sin¬ 
gle  loop  system,  digital  design  is  very  difficult,  and  especially  there  is  no  sufildait  the¬ 

ory  for  a  multi -loop  feedback  system,  the  hybrid  dectrohydrauHc  servo  system  is  widely 
used  in  almost  all  the  control  system  design  of  the  table  No  doubt,  this  scheme  is  unfavorable 
to  improving  its  accuraqr  and  rdiat^ty.  Of  caurse;  optimal  control  and  adaptive  control  are 
convenient  for  digital  implementation,  but  they  will  lake  a  long  time  to  be  carried  out  in 
software  and  this  can  not  satisfy  the  expected  bandwidth  of  a  system.  In  addition,  since  the 
wordlength  of  a  microprocessor  is  constant,  a  ooirq>lex  control  algorithm  will  make  large  num¬ 
bers  of  accumulated  errors.  In  order  to  avoid  the  above  difiiailties,  a  simple  algorithm  based 
on  a  state  observer  and  state  feedback  is  given  in  this  paper.  The  digital  simulation  results  indi¬ 
cate  that  its  performance  is  satisfactory. 

1  Algorithm  of  the  Controller 

In  order  to  obtain  acoderation  signals  and  inq)lement  the  system  convenienfly  by 
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ft  microprooessor  ,  ft  state  observer  is  used  to  gain  all  the  states  in  this  design.  The  system 
block  diagram  is  shown  in  Fig.  1. 

The  state  space  representation  of  the  system  is  ^vcn  by 

X’‘AX  +  BU 

Y^CX  (1) 

vibac  {A,B,C}  can  be  both  oontronable  and  observaMe 

Assume  the  sampling  period  to  be  To,  then  the  discrctiation  of  Bq.  ( 1 )  is  as  foDows: 

A(k  +  1) -FA(jk)  +  Gt/(fe) 

Y{k  +  l)  =CXik+l) 

where  F-e**"'  ,C«J^  e^^^df-B  and  Eq.  (2)  can  be  both  controllable  and  observabh 

Introducing  state  feedback,  we  obtain 

Uik)  =rik) -KXi,k) 

Substituting  Eq.  (3)  into  Eq.  (2),  we  have 

Xf^k+l)  =FA(fc)  +  G[r(fc)  -Ja(k)]=[F-GiqA(k)  +GKk) 

y(fe+l)  =CA(Jt+l)  (4) 

From  the  linear  system  theory,  we  know  that  the  poles  of  the  system  can  be  assigned  to 
any  expected  position  by  selecting  different  feedbadc  matrix  K,  and  tite  introduction  of  state 
feedback  will  not  destroy  the  oontrollablity  of  the  system.  Since  Eq.  (2)  is  observable^  we  can 

introduce  a  dosed  loop  state  observa  to  obtain  all  the  states .  Generally  ,  the  observer  has 
two  forms; 

^k+l)  =  =F^k)+  GUik) 

^k  +  1)  =  F5(k  +  1) +LCnk+l) -f(jk  +  l)] 

iTk+l)=Cj^k+l)  (5) 

ft^ere  A![k)  is  the  reconstructed  state  of  A(k),  and  L  is  the  gain  matrix  of  the  observer. 

(2) 

k  +  1)  =F^k) +Glf(jk) +L[y(k) -f(k>] 
nk+i)=-c]ik+i) 

where  L  is  the  same  as  the  above-mentioned. 


(2) 


(3) 


(6) 


Substituting  J^k)  into  Eq.  (3),  we  have 

Uik)  -KXik)  (7) 

Substituting  Eq.  (7)  into  Eq.(  2),  we  obtain 

A(k+1)  =FA(k) -GX^k) +GKk) 
y(k  +  l)=CA[k+l) 

In  the  continuous -time  system,  the  separation  theorem  shows  that  the  dynamics  of  the 
plant  is  indqjendait  of  the  dynamics  of  the  observer,  and  vice  versa.  For  a  discrete -time  sys¬ 
tem,  there  is  also  a  same  theorem. 

Let  AX(k)=2t(k)-ik),  then 

Xik  +  l)  ={F-GK)Xik)  +GKAX{k)  +GKk) 
AX(k+l)=(F-LC)AX(k) 

Rewrite  Eq.  (9)  in  the  vector  form  as  follows 


Xik  +  l) 

AA(k  +  l). 


(10) 


From  the  matrix  theory,  the  duuractcristic  polynomail  of  the  disetde  transfer  matrix  of 
Eq.  ( 10)  is  pven  by  Eq.  ( 11) 

det  j-(;i  I-F+GK)  (A/-  F+LC)]  =  dct[A/-F-l-GX]  det  [A/-F+LC]=0  ( H) 

So  the  poles  of  the  plant  and  the  observer  can  be  assigned  to  the  expected  position  by 
scparatdy  selecting  X  and  L.  Consequently  the  controller  and  the  observer  can  be  designed 

sqaratdy.  u  a.  i 

The  stability  of  the  system  is  discussed  as  follows.  From  Eq.  (2),  so  long  as  the  whole 

eigenvalues  of  the  transfa  matrix  are  inside  the  unit  drde  ,  Eq.  ( 10)  is  stable.  When  the  system 

comes  into  a  steady-state,  Eq.  ( 12)  is  obtained. 

lim^Wk)  =lim(F-LC)‘ AA(0)  =0 

Eq.  ( 12)  shows  when  k  is  sufficiently  larger  the  state  error  of  the  observer  and  the  plant  is 
zero,  i  e. ,  the  states  estimated  by  the  observer  are  the  real  states  of  the  plant 

Generally,  the  whole  state  initial  values  of  the  plant  arc  impossffile  to  be  known,  so  , 
when  X(0)  #  0 ,  at  the  beginning  of  the  systan  transient  procedure  the  observed  states  are 
not  the  actual  states  of  the  plant;  and  when  A(0)  =0,  the  observed  states  are  the  actual  ones 
of  the  plant  at  the  beginning.  Analysing  Eq.  (6).  we  find  that  if  X(0)  =0.  then  the  observer 
matrix  L  is  useless  and  an  undosed-loop  obsenra  will  satisfy  the  requirements. 
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2  Pure  Digital  Control  for  the  Outer  Axis  of  a  Three  ^  Axis  Flight 
Simulation  Table 

The  dectrohydiaulic  servo  system  used  in  the  outer  axis  of  a  throe~axis  table  can  be  ec- 
ptessed  approximatdy  as  follows. 

S'  1  * 

&  •  -g-  .0), 

5(sH2<*ai,s +0)  ) 

=800m]/rad  — displacement  of  the  motor 
CD,  =207.35iad/s  —  characteristic  angular  fiequen^  of  the  motor 
{=0.1— danqring  ratio  of  the  motor  - 

K —  gain  of  ampliiler  (determined  by  speed  errors) 

The  following  pioblentis  will  be  considered  in  design. 

11  Sdecting  Expected  Poles 

From  Eq.  ( 13) ,  we  know  that  flte  system  is  of  three -order.  Then,  we  can  make  flie  ex¬ 
pected  poles  of  the  jdant  equal  to  the  dosed— loop  pedes  of  a  three  -order  optimal  system. 
The  expected  system  is  expressed  as  follows. 

“  5^+2<»,sH2<u?5+eD*  ^ 

where  cD,=lif„,  and  is  the  open  loop  gain  of  the  system 

Her^  we  take  ^1,= 50  and  substitute  it  into  Eq.  ( 14) ,  and  then  we  obtain  the  expected 
poles  of  the  plant  as  follows. 

2,  =  - 106.0639,  Ay  = -47 ± /85 

In  order  to  obtain  the  actual  states  of  the  plant  as  fast  as  possible;  the  expected  poles  of 
the  observtt  should  be  as  lar  away  from  the  Im  axis  in  [s]  plane  as  possible,  but  in  jwactice 
it  Is  limited  by  the  osdllating  frequency  of  the  dodc.  Generally  its  real  parts  are  takf"  as 
(2  ~  3)  times  as  the  real  parts  of  the  expected  poles  of  the  idant  Here,  we  take 
a;=21,= -212,  .l^=2ReCA,a]±740 

2.2  Sdectii^  the  Samplii^  Paiod 

The  sampling  period  To  must  satisfy  the  san^ling  theorem.  When  is  too  'large;  some¬ 
times  we  can  obtain  an  expected  robust  discretization  system,  but  after  state  feedback  being  in¬ 
troduced  the  dosed— loop  system  will  be  an  unstable  positive  feedback  system  ( the 
undoscd-loop  system  must  be  stable  ) .  When  T,  is  too  Kttle;  the  accumulated  oror  will 
make  the  S3^tem  unstable  and  reduce  its  accura^. 

So,  from  the  viewpoint  of  engineering,  To  wOl  be  sdected  as  follows. 
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0)  •  >  niax(  5~10)  x( or  the  bandwidth  of  the  system) 

In  this  paper,  the  charactenstic  hequen^  of  the  observer  is  larger  than  that  of  the  ex¬ 
pected  system,  and  it  has  decisive  influence  on  the  system  So  to,  should  be  the  diaractcristic 
iiequency  of  the  observer.  ITien 


Here,  Tq  is  sdected  as  O.OOOSs. 

2.3  CoDStraction  of  the  Control  3ystan 

The  inplementation  block  diagram  is  shown  in  Fig.  2  The  controller  and  the  observer  in¬ 
side  the  dosed  dotted  line  are  iaq>lemented  by  a  microprocessor.  Because  the  introduction  of 
state  feedback  dan^  the  output  of  the  system,  we  introduce  a  conqsensator  into  the  input 
pass.  Its  parameters  are  determined  as  follows: 

From  the  viewpoint  of  the  continuous  system  and  assuming  that  there  is  no  obsover,  the 
system  shown  in  Fig.  2  is  equivalent  to  the  system  shown  in  Fig.  3.  Because  ik,  is  the 
gain  of  state  x,  and  considering  that  F-c  i  x ,  ,flic  ou^ut  feedback  gain  in  Fig.  3  should  be 
k,/c, .  For  the  unit  step  input,  the  ouQiut  of  the  system  is 

y(5)=- - ^ - .4-  = - ^ - - 

(sG'(s)  +  -^c,)  ^  sG(s)+kt's 

c  I 

Then  y(«®  )  “Km  s.  l^s)  =c,/fc, 

•—0 

So,  in  order  to  ensure  y(oo  )  =1,  the  block  ?  in  Fig.  3  should  be  rqilaoed  by  an  amplifer 
with  a  gain  Jlc,/c ,. 

3  Simulation  Results 

Based  on  Bq.  ( 13) ,  digital  simulation  has  been  done  The  stqi  and  sine  response  of  the 
frfant  with  rero  and  little  and  large  initial  states  are  discussed  lespectivdy.  When  the  state  ini¬ 
tial  values  of  the  plant  are  zero,  the  system  is  equal  to  one  without  an  observer.  Its  response  is 
shown  in  Fig.  4.  The  unit  stq>  response  of  the  {dant  with  initial  states  is  shown  in  Fig.  5. 
When  the  state  initial  values  of  the  plant  is  not  zero,  a  non -normal  balance  point  occurs  in 
tiie  position  following  time.  Conqiaring  Fig.  4  with  Fig.  5,  we  find  that  the  position  following 
time  is  about  O.OSs  and  0.1s  respectivdy  for  the  system  without  and  with  an  observer.  The  sin 
(79t)  response  (for  the  gain  fiequengr  bandwidth  test)  and  the  sin  (69t)  response  ( for  the  phase 
frequency  banduidth  test )  of  the  system  are  shown  in  Fig .  6a  and  Fig .  fib  respectivdy. 

4  Conclusions 

The  pure  digital  control  system  for  the  outer  axis  of  a  three— axis  flight  simulation  table  is 
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designed  in  this  psper.  It  is  single  and  convenient  for  implementatinn,  Simulation  results  show 
that  although  a  non— normal  balance  point  occurs  in  the  system  with  an  observer,  the  jerhmw 
given  in  this  paper  can  ensure  perfect  dynamic  performance  for  the  three— axis  table  and  suit 
for  a  most  dectrohydraulic  position  servo  systoo. 


Fig .  3  Ihe  equivalem  system 
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Rg .  4  The  unit  step  response  of  the  plant 
with  a  2eio  initial  value 


(a) 


0  ao5  at  ai5  02  a25 


Fig .  S  The  unit  step  response  of  the  plant 
with  a  little  in^  value 


1— plant  2— ofisaver  3 — input 


- - - I 

0  ao5  at  ai5  a2  a2s 


TIME(SEC) 

(b) 


Fig .  6  The  sinewave  response  of  the  {dant  with  an  initia]  value 
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COMPUTERS 


A  Hybrid  Learning  Method  for  Multilayer 

Perceptrons 

OF  HARBIN  INSTITUTE  OF 

ui?  ifctHNOLOGYJ  in  Zhao  Made  Huang  Wenhu  mhgjiarong 

(School  oT  Astronautics) 


Alistract  A  Newton  learning  method  for  a  neural  network  of  multilayer  perceptrons  is  pro¬ 
posed  in  this  paper.  Furthermore,  a  hybrid  learning  method  is  legitimately  devdoped  in 
combination  of  die  backpropagation  method  proposed  by  Ruirelfaart  et  al  with  die 
Newton  learning  method.  Finally,  the  hybrid  learning  algorithm  is  compared  with  the 
backpropagation  algorithm  by  some  illustrations,  and  the  results  show  that  this  hybrid 
learning  algorithm  has  the  characteristics  of  rapid  convergence 
Key  Words;  Computer;  Newton  learning  method;  naval  network;  multilayer  perception 


0  Introduction 


Multilayer  peioqitiDns  are  feed -forward  neits  with  one  or  more  layers  of  nodes  between 
the  input  and  output  nodes.  These  additional  layers  oontairi  hidden  units  or  nodes  that  are 
not  directly  coimected  to  both  the  input  and  output  nodes.  Multilayer  perceptrons  overcome 
many  of  the  limitations  of  single  layo’  pero^trons  ^  ,  but  were  generally  not  used  in  the 
past  because  effective  training  algorithms  were  not  available.  This  has  changed  with  the  devel¬ 
opment  of  new  training  algorithm,  L  e.,  the  backpropagation  algorithmic.  Since  its 
emergmoe,  they  have  been  quite  successfully  used  for  a  wide  variety  of  tasks  such  as  informa¬ 
tion  processing,  coding  and  pattern  recognition  .  Unfortunatdy,  in  many  applications, 
the  number  of  iteration  required  before  convergence  can  be  large  and  there  exists  local  minima 
problem.  Some  mcxlifications  to  the  backpropagation  algorithm  were  given  by  Yogi  et  al^  in 
order  to  accelerate  the  convergence  of  the  backpropagation  algorithriL  But  we  found  that 
these  modifications  caruiot  accelerate  the  convergence  by  many  numerical  simulations,  hi  Sec¬ 
tion  2  of  the  paper  the  backpropagation  algorithm  and  some  modifications  to  it  are  briefly  des- 
’  cribed.  In  Section  3  tiie  recursive  formulas  for  calculating  Hessian  matrix  of  multilayer 
peroqitrons  are  derived  and  a  Newton  learning  algorithm  is  given.  A  hybrid  learning  algorithm 
for  multilayer  percqitrons  is  proposed  in  Section  4,  and  'in  Section  5  two  examples  are 
given  to  oonqiaie  tile  two  methods. 

1  The  Backpropagation  Method  and  Its  Modification 

LI  Modd  of  Neural  Nets 


ito 


In  general,  a  neural  network  consists  of  N  nodes  (processors),  each  of  which  is  con¬ 
nected  to  all  the  other  nodes.  In  the  present  work,  however,  only  layered  feedforward  networks 
are  used  .Use  of  such  a  layered  structure  sacrifices  a  certain  amount  of  parallelism  since  only 
processors  in  the  same  lay®  operate  in  parallel  but  adjacent  layers  compute  senally.  Although 
some  parallelism  is  lost,  these  networks  are  popular  because  there  exists  a  sophisticated  learn¬ 
ing  scheme  for  training  the  network.  A  typical  network  using  tl’e  layered  structure  is  shown  in. 
Fig.  1.  It  consists  of  an  input  layer,  an  output  layer,  and  hidden  layers  in  between.  The  nodes 
in  the  hidden  layer  arc  necessary  to  implement  nonlinear  mappings  between  the  input  and 
output  patterns. 

Associated  with  each  node  is  a  threshold  value.  The  output  of  each  node  is  determined 
by  the  levd  of  the  input  signal  in  relation  to  the  threshold  value.  Signals  can  propagate 
through  the  network  only  between  adjacent  layers  and  only  in  a  single  direction.  These  sig¬ 
nals  are  naodified  by  the  connection  strengths  (.weights)  between  the  ncxles.  Let 
threshold  of  the  tth  node  in  the  kth  layer, 

IFj  j  t  weight  betweenjth  node  in  the  (k  —  l)th  layer  to  ith  node  in  kth  layer, 
nett  ,k  “Piit  to  file  rth  node  in  the  kth  layff, 

Oj  n  output  of  the  tth  ncxle  in  the  kth  layer. 

The  input  to  a  ncxie  is  given  by 


The  summation  index  J  extends  over  all  nodes  in  the  ( k  —  l)th  layer.  The  output  of  a  node  is 
given  as 

=;(  net  iO  (2) 

where  f  is  a  sigmoid  function  given  by 


1.2  The  Backpropagation  Method  and  Its  Modification 

In  order  to  describe  conveniently,  here  we  use  vector  notation .  The  basic  iteration  formu¬ 
la  of  the  backpropagation  algorithm  is  : 

where  IF  is  a  column  vector  whidi  involves  all  the  weights  and  the  thresholds,  W*  rqjresents 
the  column  vector  of  IF  in  kth  iteration,  £(  IF)  is  the  total  error,  I  e. ,  the  sum  of  the  squares 
r-f  the  error  of  the  ouq)ut  signal ,  V  £(  IF)  is  the  gradient  colu  iin  vector  of  £(  IF)  with  respect 
to  IF,  is  a  learning  constant ,  and  is  a  momentum  factor. 

In  the  backpropagation  algorithm,  rik  and  a*  are  preset  constants  .  Let  £  ( IF)  =JE^  IF), 
where  £p(  IF)  is  the  error  function  under  the  pth  input  patterns,  I  c. , 


where  the  summation  index  j  extends  ovct  all  nodes  in  the  ou^ut  layer,  dpj  is  the  desired  out¬ 
put  ( provided  by  a  teacher  via  the  input -output  samples)  at  the  ith  output  node  and  is 
the  actual  ouQjut  at  the  jth  output  node  The  element  of  V  E(  W)  is : 

^Aik  ®  for  ““f  i  of  the  Irth  layer  under  the  pth  input  samples. 

J  prk=n 

i  Z(5py.it+,  k+i)  otherwise 

wiicrc  the  index  s  indicates  summation  over  all  nodes  r/  the  kth  layer.  The  detailed  desoiption 
of  the  backpropagation  algorithm  is  given  in  [2]  and  ;^4j. 

Some  modifications  to  the  backpropagatioii  algorithm  was  given  in  [5].  These  modifica¬ 
tions  are  as  follows:  Instead  of  updating  the  network  weights  after  eadi  pattern  is  presented 
tc  the  network,  the  network  is  updated  only  after  the  entire  repertoire  of  patterns  to  be  learn¬ 
ed  has  been  presented  to  the  network;  Instead  of  keeping  tj  constant,  it  is  varied 
dynamically  so  that  the  algorithm  utilizes  a  near -optimum  n  ;  the  momentum  factor  a  is 
set  to  zero  when,  an  signified  by  a  failurc  of  a  step  to  reduce  the  total  error,  the  information 
inherent  in  prior  stqis  is  more  likdy  to  be  misleading  than  beneficial.  Only  after  the  network 
takes  a  useful  step,  i  e. ,  one  that  reduces  the  total  error,  does  a  again  assume  a  non— zero 
value. 

Here,  we  propose  a  new  scheme  for  selecting  x^  and  they  are  selected  sudi  that 
(^k,Vk  =  nrnEiW^-nV  EiW'‘)+oi[W'‘-W'‘-’)) 

Theoretically  we  have  solved  the  problem  of  how  to  select  a  and  i,  in  the  backpropagation 
algorithm.  However,  by  many  numerical  simulations  we  find  that  both  the  modrications  to 
the  badepropagation  algorithm  given  in  ' [5]  and  above  sdieme  cannot  nnyiirratf  the 
convergence  of  the  backpropagation  method.  In  the  next  section  we  propose  a  new  learning 
method  for  multilayer  pero^trons. 

2  The  Newton  Learning  Method  for  Multilayer  Perceptions 

In  order  to  solve  above  problems,  we  propose  a  Newton  learning  method  for 
peroeptrons. 

2.1  The  Hessian  Matrix  of  M»tltiiayer  Peroeptrons 

The  Hessian  matrix  of  the  network  under  the  pth  input  patterns  is  denoted  as 
we  have 


multilayer 


Hp.  Hence 


(5) 

(6a) 

(6b) 


where  V*£,  is  a  second  orda  derivative  matrix  of  £,  with  respect  to  W.  For  convenience;  we  in¬ 
troduce  the  following  notation: 


H'”  =  — 

L  1  I  ^ 


ffl  » 

a««,. 

Since  H,  is  a  symmetric  matrix,  we  can  discuss  it  in  two  cases : 
l)Thc  diagonal  dements  of  H, 

From  (5),  we  have; 

. =-5  0 

By  (7a),  we  obtain 


!FE, 


tr _ r_:a _  n  dd^i. 


smoe 


-P^sftr..  =  Snetfj.f 

OWjj.r  ^n.r 


(7a)  • 


(7b) 

(7c) 

(7d) 


According  to  (1)  and  (7b),  we  can  obtain 


SWj,l,r 


Hence 


~  yij.r 


2) The  non -diagonal  danents  ofH,  i  (r  >q) 
From  ( 7a)  and  ( 5) ,  we  have 

If  '’P  _  ^  Snet , 

anet;.,.,  ’  awj.,.. 


_  8net,,„  , ^  d5^.,  a c  . 


(Sa) 


According  to  ( 1) ,  (7c)and(7d),  we  can  obtain 


w4icre  the  lecuisive  foimulas  for  calculating  and  j.,,,  are  given  in  Appendix. 


2.2  TTie  Newton  Learning  Method  for  Maltilayer  Peroeptrons 

According  to  the  Hessian  matrix  derived  above,  we  can  obtain  the  basic  iteration  formula 
of  the  Newton  learning  method  for  multilayer  peroq>trons  as  follows: 


^  G,  (9) 

where 

r  p 

V£;(W*) 

p 

The  index  p  represents  summation  over  all  input— output  samples  used  for  training,  ^  is  a 
constant  whidi  is  called  a  Newton  lea  ning  constant  in  the  paper. 

The  Newton  learning  method  has  been  implemented  on  a  Sun3  Workstation  in  Clanguage. 
We  find  by  many  illustrations  that  this  algorithm  cannot  converge  to  a  solution.  There  are  pos¬ 
sibly  severd  reasons  for  it  through  our  analysis:  1 )  Although  the  Newton  algorithm  can 
quadratically  converge  near  the  exact  solutions,  it  is  very  difficult  to  find  a  better  initial  guess; 
2)Aiter  several  iterations,  the  matrix  H  will  become* singular;3)We  cannot  ensure  that  for  each 
iteration  the  quadratic  model  can  adequatdy  approximate  the  objective  function  E. 

To  apply  the  Newton  learning  method  proposed  in  the  paper,  we  propose  a  hybrid  learn¬ 
ing  method  which  combines  the  backpropagation  method  with  the  Newton  learning  mdhod. 
The  hybrid  learning  method  has  the  advantages  of  the  two  methods.  The  illustrations  given  in 
Section  5  show  that  this  new  learning  algorithm  has  a  surprising  speed  of  convergence. 

3  The  Hybrid  Learning  Method  for  Multilayer  Perceptrons 

Since  the  backpropagation  method  is  a  steepest  decent  mdhod,  it  can  globally  converge  to 
a  solution.  The  Newton  learning  method  can  quadratically  converge  in  a  local  region  of 
the  solution.  Therefore,  the  hybrid  learning  method  can  rapidly  globally  converge.  The 
basic  idea  of  the  hybrid  learning  method  is  :  after  each  input  pattern  is  presented  to  a  net¬ 
work,  the  weights  and  thresholds  are  updaled-^by  the  backpropagation  algorithm;  after  entire 
input  patterns  are  prcsaitcd  to  the  network,  the  weights  and  thresholds  are  updated  by  the 
Newton  learning  method  The  detailed  description  of  the  algorithm  is  as  follows. 

Stepl  Initialize  weights  and  thresholds  ( set  all  weights  and  node  thresholds  to  small  ran¬ 
dom  values). 

Stcp2  According  to  (1)  and  (2),  calculate  the  actual  outputs  of  each  nede  under  ui: 
pth  input  pattern  in  forward  direction. 
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Stcp3  Calculate  the  error  signals  of  the  output  nodes  in  terms  of  ( 6a)  and  update  the 
weights  and  thresholds  of  the  output  layer  aocoiding  to  ( 3) 

Step4  Calculate  the  ,  of  the  ou^ut  layer  in  terms  of  (  A.  la)  and  calculate  the  corre¬ 
sponding  elements  of  the  matrix  by  ( 8a) . 

Steps  Calculate  and<p,,y^f^,m  terms  of  (A-3)  and  (.A 4),  then  calculate  the  corre¬ 
sponding  elements  of  the  matrix  by  ( 8b) . 

Stq)6  By(Al  )calculate  the  y,  ^ ,  of  the  penultimate  nodes,  and  calculate  the  corresponding 
dements  of  the  matrix  ff,  by  ( 8a) . 

Step?  Calculate  the  error  signals  of  the  penultimate  nodes  in  terms  of  ( 6b) ,  and  update 
the  weights  and  thresholds  of  the' penultimate  layer  aax)rding  to  ( 3) . 

Step8  Repeat  Step  S,  Step  6  and  Step  7,  until  the  entire  elements  of  matrix  is  calcu¬ 
lated. 

Step9  Repeat  Stqp  2— Step  8,  until  all  input  patterns  are  presented  to  the  network. 

.  SteplO  Decide  uiiethcr  or  not  matrix  H  has  an  inverse  matrix  .If  it  has,  then  go  to  Step 
11,  otherwise  go  to  Step  12 

Stepll  Update  all  weights  and  thresholds  of  the  network  according  to  ( 9) . 

Stepl2  If  the  maximum  of  the  dilTetences  between  the  actual  output  and  the  desired  out¬ 
put  of  the  output  layer  for  all  input  patterns  is  less  than  preset  value,  then  termi- 
nate^  otherwise  go  to  Step  2 

4  Illustrations  and  Conparison 

The  backpropagation  algorithm  and  the  hybrid  algorithm  have  been  implemented  on 
Sun3  Workstation  in  C  and  FORTRAN  languages.  Two  examples  are  given  as  follows. 

4.1  The  XOR  Problem 

The  XOR  problem  is  a  dilHcult  dassiral  problem,  which  can’t  be  solved  by  single  layer 
peroeptrons.  Since  it  is  a  simplest  nonlinear  seperatabte  problem.  The  XOR  problem,  L  e  ,  Ex- 


dusive  OR  problem  is  the  problem  that  a  neural  net  can  generate  the  following  outputs  for  the 
following  input  samples  by  training. 

Input 

Pattern 

Output 

0 

0 

0 

1 

1 

0 

0 

1 

1 

1 

0 

1 

Fig.  2  and  Fig.  3  show  the  solution  solved  by  the  badepropagation  algorithm  and  the  hybrid 
algorithm  respectivdy.  The  solution  in  Fig.  2  is  reached  afler  532  sweeps  through  the  four 
stimulus  patterns  with  learning  rate  of  ij  =0.5.  and  momentum  factor  of  a  =0.9. The  solution 
in  Fig.  3  is  reached  after  39  sweeps  with  tj  =0.5,  a  =0.9  and  tj/  =0.11. 

Through  many  numerical  simulations  we  (Ind  the  parameter  rt  and  a  effects  bn  the  num¬ 
ber  of  iterations  before  convergence  for  the  two  methods.  Fig  .4  and  Fig.  5  show  these  effects 


respectlvdy.  The  Newton  learning  constant  \l/  affects  the  number  of  iterations  before 
convergence  for  the  hybrid  algorithm.  Fig  .6  shows  this  effect  ip  is  now  varied  between  0.005 
and  0.02  by  ocperience. 

In  numerical  simulation  we  also  find  that  the  Newton  updation  for  the  weights  and 
thresholds  only  occur  in  the  first  several  sweeps  in  the  hybrid  algorithm  Table  1  is  ob¬ 
tained  by  assuming  that  after  10  sweeps  the  Newton  updating  stops. 


Table  L  Set  a=0.95and  |/f=0.011 


n 

Backpropagation 

Hybrid  Learning 

sweep 

CPU  time 

sweep 

CPU  time 

0.1 

2736 

246 

127 

6 

0.2 

1011 

87 

124 

6 

0.3 

581 

51 

68 

4 

0 

383 

35 

29 

2 

0.5 

274 

25 

^ . . 

28 

2 

BP — Badcpcopagation  HL — Hy  bno  Learning 


4.2  TTie  Parity  Pioblan 

The  parity  problem  is  a  problem  iriiidi  can  computer  the  parity  of  any  input  pattern,  i 
t ,  the  output  required  is  1  if  the  input  pattern  contains  an  odd  number  of  Is  and  0  other¬ 
wise  This  is  a  very  difficult  problem  because  the  most  similar  patterns  require  different  an¬ 
swers.  The  XOR  problem  is  a  parity  problem  with  input  patterns  of  size  two.  Here  we  have 
done  a  parity  problem  with  patterns  of  size  3.  Fig.  7  and  Fig.  8  show  the  actual  solutions  with 
three  input  lines  and  three  hidden  units  obtained  by  the  backpropagation  algorithm  arid  the  hy¬ 
brid  algorithm  respectively.  Let  n  =0.5,  a  =6.9,  and  tj,  =0.006.  The  solutions  in  Fig.  7  and 
Fig .  8  are  reached  after  1441  and  52  sweeps  through  the  eight  stimulus  patterns  respectively. 

5  Conclusions 

In  this  paper,  a  Newton  learning  algorithm  for  multilayer  peiceptrons  is  derived  and  a  hy¬ 
brid  learning  algorithm  is  devdoped  in  combination  of  the  backpropagation  algorithm  with 
the  Newton  learning  algorithm  .  The  results  of  numerical  simulations  given  in  the  pape.  shov% 
that  the  hybrid  algorithm  has  much  faster  convergence  than  the  backpropagation  algoritima 
The  future  study  of  the  project  is  to  generalize  this  new'  algorithm  to  high  order  ncuia,  noLsand 
apply  it  to  solving  practical  problems. 
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Appendix 

TTic  Rcairsive  Calculation  For  and 

1.  The  Recursive  Calculation  For  7^^, 

We  derive  the  fonnula  of  yj,  in  two  cases  as  follows. 

( i)  Layer  r  is  an  ou^ut  layer 
From  ( 6a) ,  we  have 


^pi,r  ) 

By  (7b),  we  can  obtain. 


\j  Xnet,t,,)idfi-Oft  ,)  -[(y '(nfit,i.,)y .  othayvise 

( ii)  Layer  r  is  a  hidden  layer 
From  (6b) ,  we  have 


^pt,r  J pt.r)^  ^pj  ,ri-l 

J 


Wj 


j  pi,r)  S  +I  * 

1 


^Pj.r 


By  (7c) ,  we  obtain. 


(Ala) 


Jlj.r  y  J 

(2)Fori=j. 

V.J.. -J  S„,„  w . .,)  +  j  f. , 

J  s  *  ‘  * 

2 .  The  Recursive  Cafoilation  For  Ci,j„  ( Ale) 

According  to  ( 7c) ,  we  have . 


dS, 


dnet 


r, ,  =  -y^pf.r  =y/ 

^Pj.i  .  Snet,,^^,  dnet,j,, 

..  _JOpi^  ^  r,  ^  ^ 

and  V  -  =  I 

1  r 


(A2) 


,  jor  r=q  +  l 
,jor  r>q  +  l 


_f  jbrr=^+J 

••  ti.j.r.q  -  -j  • 

J  ,  pr  r>q+l 

3.  The  Reairsivc  Calculation  For  o  -  • 

(i)Forr=q, 

By(7d),  wc obtain. 


SO„i  r  ^\0,pri9ij 
^^n.f  XiXnet^j  ,)  pr  i^j 


(u)Forr  >  q 

According  to  ( 7d)  and  by  diain  rule ,  we  have 


=  —'T'  OWgty,.^ 

^^pj.f  t  dnetft^f^;  dnet 


By  (A.  2) ,  we  obtain 


<Pi.J.r.,  =  7  t  )  -I  ( 


Output  layer 


layer  (^3^  CZ^ 


uth  nodO 


c>i 


(A:-l)th  layer ^  ^ 

V  V 


Hidden  laver 


Input  layer 


Fig .  1  The  stnicture  of  'a.  fecd-forwaid  net 
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